Ursus Americanus - Elevation Analysis

# Loads data
load("bc_bear_occurrences.Rda")
load("BC_Covariates.Rda")

# Generates summary of loaded data
summary(DATA)

##            Length Class           Mode
## Window      1     SpatialPolygons S4  
## Elevation  10     im              list
## Forest     10     im              list
## HFI        10     im              list
## Dist_Water 10     im              list

Visualizes bear occurrences in BC

# Extracts location columns
bears_loc <- occ_data[, c("decimalLongitude", "decimalLatitude", "month", "year")]
bears_loc_filtered <- subset(bears_loc, year %in% c(2020, 2021, 2022, 2023, 2024))

# Creates sf object
bears_sf <- st_as_sf(
  bears_loc_filtered,
  coords = c("decimalLongitude", "decimalLatitude"),
  crs = 4326  # WGS84 (longitude/latitude)
)

# BC Albers projection
bears_sf_proj <- st_transform(bears_sf, crs = 3005)

# Extracts coordinates
coords <- st_coordinates(bears_sf_proj)

# Extracts BC window
window_sf <- st_as_sf(DATA$Window) 
window_proj <- st_transform(window_sf, crs = 3005) 
window <- as.owin(window_proj) 

# Creates ppp object
bears_ppp <- ppp(
  x = coords[,1],
  y = coords[,2],
  window = window
)

## Warning: 187 points were rejected as lying outside the specified window

## Warning: data contain duplicated points

# Plots bear occurrences
plot(bears_ppp, 
     pch = 21, 
     main = "Black bear occurrences in BC, 2024")

## Warning in plot.ppp(bears_ppp, pch = 21, main = "Black bear occurrences in BC,
## 2024"): 187 illegal points also plotted

Compares seasonal distributions

# Defines the seasons and corresponding colors
seasons <- list(
  winter = c(12, 1, 2),
  spring = c(3, 4, 5),
  summer = c(6, 7, 8),
  autumn = c(9, 10, 11)
)

# Creates an empty list to store the ppp objects for each season
ppp_list <- list()

# Creates an empty list to store the filtered data for each season
bears_sf_list <- list()

# Loops over seasons
for (i in 1:length(seasons)) {
  
  # Filters for each season
  season_name <- names(seasons)[i]
  season_months <- seasons[[i]]
  
  # Filters bears_loc_filtered by the season's months
  bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% season_months, ]
  
  # Prints the number of observations for this season
  cat(season_name, ": ", nrow(bears_loc_season), " observations\n", sep="")
  
  # Creates sf object for the filtered season
  bears_sf_season <- st_as_sf(
    bears_loc_season,
    coords = c("decimalLongitude", "decimalLatitude"),
    crs = 4326  # WGS84 (longitude/latitude)
  )
  
  # Transforms to BC Albers projection
  bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
  
  # Stores sf object in the list
  bears_sf_list[[season_name]] <- bears_sf_season_proj
  
  # Extracts x, y coordinates
  coords <- st_coordinates(bears_sf_season_proj)
  
  # Creates ppp object for each season
  suppressWarnings(
    bears_ppp_season <- ppp(
      x = coords[, 1],
      y = coords[, 2],
      window = window
    )
  )
  
  # Stores ppp object in the list
  ppp_list[[season_name]] <- bears_ppp_season
}

## winter: 76 observations
## spring: 808 observations
## summer: 1862 observations
## autumn: 841 observations

# Plots four maps in a 2x2 grid
par(mfrow=c(2,2), mar=c(1,1,1,1)) 

# Loops over ppp objects and plots
for (i in 1:length(ppp_list)) {
  season_name <- names(ppp_list)[i]
  suppressWarnings(
    plot(ppp_list[[season_name]], 
         main = paste(season_name, "(2024)"), 
         pch = 21)
  )
}

# Visualizes covariates
plot(DATA$Elevation)

plot(DATA$Forest)

plot(DATA$HFI)

plot(DATA$Dist_Water)

# Isolates covariates
elev <- DATA$Elevation
cover <- DATA$Forest
dist_water <- DATA$Dist_Water
hfi <- DATA$HFI

Estimates the intensity of bear occurrences based on elevation

# Intensity as a function of elevation overall
rho <- rhohat(bears_ppp, elev) 
plot(rho, 
     xlim=c(0, max(elev)), 
     main="Estimated Rho vs. Elevation")

# Intensity as a function of elevation broken up by seasons
par(mfrow=c(2,2), mar=c(2,2,2,2))

for (i in 1:length(ppp_list)) {
  season_name <- names(ppp_list)[i]
  rho <- rhohat(ppp_list[[season_name]], elev) 
  plot(rho, 
       xlim=c(0, max(elev)), 
       main = paste(season_name))
}

# Loads elevation raster
elev <- DATA$Elevation

# Converts im to raster
im2raster <- function(im_obj) {
  m <- as.matrix(im_obj)
  r <- raster(m)
  ext <- c(im_obj$xrange[1] - im_obj$xstep/2,
           im_obj$xrange[2] + im_obj$xstep/2,
           im_obj$yrange[1] - im_obj$ystep/2,
           im_obj$yrange[2] + im_obj$ystep/2)
  extent(r) <- ext
  return(r)
}
elev_raster <- im2raster(elev)

# Extracts elevation values at occurrence locations
bears_sf_proj$elev_value <- extract(elev_raster, bears_sf_proj)
summary(bears_sf_proj$elev_value)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   11.68  461.09  644.88  754.84  937.91 2308.18    1499

KDE to estimate the concentration of bear occurrences based on elevation values. The elevation data for bear locations ranges from a minimum of 11.68 meters to a maximum of 2308.18 meters, with most occurrences concentrated between 461.09 meters and 937.91 meters. The mean elevation is 754.84 meters, slightly higher than the median, indicating a skew towards higher elevations.

# KDE for bear locations
kde_b <- density(bears_sf_proj$elev_value, na.rm = TRUE, bw = "SJ-dpi")

# Summarizes
summary(bears_sf_proj$elev_value)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   11.68  461.09  644.88  754.84  937.91 2308.18    1499

plot(kde_b,
     main = "Kernel Density Estimate of Elevation at Bear Locations",
     xlab = "Elevation (meters)",
     ylab = "Density",
     col = "#2C6B2F", 
     lwd = 2)  

Models elevation

# NA with mean elevation
mean_elev <- mean(as.vector(as.matrix(elev)), na.rm = TRUE)
elev_clean <- eval.im(ifelse(is.na(elev), mean_elev, elev))

# Establishes elevation model
elevation_model <- ppm(bears_ppp, ~ elev, covariates = list(elev = elev_clean))
print(elevation_model)

## Nonstationary Poisson process
## Fitted to point pattern dataset 'bears_ppp'
## 
## Log intensity:  ~elev
## 
## Fitted trend coefficients:
##  (Intercept)         elev 
## -17.36871576  -0.00253243 
## 
##                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -17.36871576 2.963479e-02 -17.426798885 -17.310632643   ***
## elev         -0.00253243 4.234993e-05  -0.002615434  -0.002449426   ***
##                   Zval
## (Intercept) -586.09208
## elev         -59.79774

# Predicts and plots intensity
pred_intensity_elev <- predict(elevation_model)
plot(pred_intensity_elev, 
     main = "Predicted Intensity Based on Elevation")
points(bears_ppp, 
       pch = 20, 
       cex = 0.4,
       col = "#B24A2F")

Fitted model: Intercept: -17.37; Elevation coefficient: -0.0025. The negative coefficient for elevation indicates that, as elevation increases, the expected intensity of bear occurrences decreases. This suggests that bears are less likely to be found at higher elevations.

Computes and plots residuals of the point process model fitted to elevation data. Using Pearson residuals and removing any non-finite residuals.

# Establishes elevation residuals
elevation_residuals <- residuals(elevation_model, 
                                 type = "pearson")

# Removes any non-finite residuals
elevation_residuals$v[!is.finite(elevation_residuals$v)] <- NA

# Plots elevation residuals
plot(elevation_residuals, 
     main = "Residuals - Elevation Model", 
     na.col = "transparent")

Partial Residuals

# Elevation model
elevation_model <- ppm(bears_ppp, ~ elev, covariates = list(elev = elev_clean))

# Computes partial residuals for elevation
par_res_elev_all <- parres(elevation_model, "elev")

# Plots partial residuals
plot(par_res_elev_all,
     legend = FALSE,
     lwd = 2,
     main = "Partial Residuals Elevation - All Seasons",
     xlab = "Elevation (km)",
     ylab = "Partial Residuals")

models_elev_seasonal <- list()

# Loops over season in the ppp_list and models elevation
for (season in names(ppp_list)) {
  models_elev_seasonal[[season]] <- ppm(ppp_list[[season]], ~ elev, covariates = list(elev = elev_clean))
}

## Warning: glm.fit: algorithm did not converge

# 2x2 plot layout
par(mfrow = c(2, 2))

# Generates and plots partial residuals per season
for (season in names(models_elev_seasonal)) {
  model <- models_elev_seasonal[[season]]
  parres_season <- parres(model, "elev")
  plot(parres_season,
       legend = FALSE,
       lwd = 2,
       main = paste("Partial Residuals - Elevation", season),
       xlab = "Elevation (km)",
       ylab = "Partial Residuals")
}

# Seasonal elevation models
models_elev_seasonal <- list()
for (season in names(ppp_list)) {
  models_elev_seasonal[[season]] <- ppm(ppp_list[[season]], ~ elev, covariates = list(elev = elev_clean))
  cat("Model for", season, ":\n")
  print(models_elev_seasonal[[season]])
}

## Warning: glm.fit: algorithm did not converge

## Model for winter :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
## 
## Log intensity:  ~elev
## 
## Fitted trend coefficients:
##   (Intercept)          elev 
## -20.281500393  -0.004896402 
## 
##                  Estimate         S.E.       CI95.lo      CI95.hi Ztest
## (Intercept) -20.281500393 0.1590970890 -20.593324958 -19.96967583   ***
## elev         -0.004896402 0.0004162485  -0.005712234  -0.00408057   ***
##                   Zval
## (Intercept) -127.47876
## elev         -11.76317
## *** Fitting algorithm for 'glm' did not converge ***
## Model for spring :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
## 
## Log intensity:  ~elev
## 
## Fitted trend coefficients:
##   (Intercept)          elev 
## -18.938423756  -0.002390713 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -18.938423756 6.296883e-02 -19.061840399 -18.815007114   ***
## elev         -0.002390713 8.708531e-05  -0.002561397  -0.002220029   ***
##                   Zval
## (Intercept) -300.75870
## elev         -27.45254
## Model for summer :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
## 
## Log intensity:  ~elev
## 
## Fitted trend coefficients:
##   (Intercept)          elev 
## -18.213839873  -0.002227723 
## 
##                  Estimate         S.E.       CI95.lo      CI95.hi Ztest
## (Intercept) -18.213839873 4.265390e-02 -18.297439975 -18.13023977   ***
## elev         -0.002227723 5.699771e-05  -0.002339437  -0.00211601   ***
##                   Zval
## (Intercept) -427.01468
## elev         -39.08443
## Model for autumn :
## Nonstationary Poisson process
## Fitted to point pattern dataset 'ppp_list[[season]]'
## 
## Log intensity:  ~elev
## 
## Fitted trend coefficients:
##   (Intercept)          elev 
## -18.514156648  -0.003139775 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -18.514156648 5.671064e-02 -18.625307457 -18.403005838   ***
## elev         -0.003139775 9.455484e-05  -0.003325099  -0.002954451   ***
##                   Zval
## (Intercept) -326.46708
## elev         -33.20586

# Combines seasonal data
bears_all_elevation <- do.call(rbind, lapply(names(bears_sf_list), function(season) {
  sf_season <- bears_sf_list[[season]]
  coords <- st_coordinates(sf_season)
  elev_vals <- extract(elev_raster, sf_season)
  data.frame(x = coords[, 1], y = coords[, 2],
             elev_value = elev_vals, season = season)
}))
bears_all_elevation$season <- as.factor(bears_all_elevation$season)
bears_all_elevation$count <- 1

Winter: Intercept: -20.28; Elevation Coefficient: -0.0049. For winter, the negative coefficient for elevation indicates that bear occurrences decrease as elevation increases. The model did not converge, so caution is needed in interpreting the results.

Spring: Intercept: -18.94; Elevation Coefficient: -0.0024. In spring, the relationship between elevation and bear occurrences is also negative, with a slight decrease in occurrences as elevation increases. However, the model showed some issues with convergence.

Summer: Intercept: -18.21; Elevation Coefficient: -0.0022. Summer shows a similar negative relationship, with bear occurrences decreasing slightly as elevation increases.

Autumn: Intercept: -18.51; Elevation Coefficient: -0.0031. In autumn, bear occurrences also decrease with elevation, with a slightly stronger negative coefficient compared to spring and summer.

All seasons show a negative relationship between bear occurrences and elevation. Winter has the largest negative coefficient (-0.0049), meaning a stronger decrease in occurrences with elevation.

# GAM with elevation and season
gam_elevation <- gam(count ~ s(elev_value) + season, data = bears_all_elevation, family = poisson())
summary(gam_elevation)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## count ~ s(elev_value) + season
## 
## Parametric coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  6.729e-17  4.510e-02       0        1
## seasonspring 2.172e-22  6.531e-02       0        1
## seasonsummer 7.375e-22  5.431e-02       0        1
## seasonwinter 8.153e-21  1.544e-01       0        1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value
## s(elev_value)   1  1.001      0       1
## 
## R-sq.(adj) =    NaN   Deviance explained =  NaN%
## UBRE = -0.99521  Scale est. = 1         n = 2088

GAM results suggest that neither elevation nor season has a significant influence on the spatial distribution of bear locations. Warnings about the model fitting process suggest that the results may not be reliable.

# Polynomial ppm model with season
bears_all_elevation <- ppp(
  x = bears_all_elevation$x,
  y = bears_all_elevation$y,
  window = window,
  marks = bears_all_elevation$season
)

## Warning: 187 points were rejected as lying outside the specified window

## Warning: data contain duplicated points

marks(bears_all_elevation) <- data.frame(season = marks(bears_all_elevation))

model_poly_elevation <- ppm(bears_all_elevation, ~ polynom(elev, 2) + marks,
                     covariates = list(elev = elev_clean))
print(model_poly_elevation)

## Nonstationary multitype Poisson process
## Fitted to point pattern dataset 'bears_all_elevation'
## 
## Possible marks: 'autumn', 'spring', 'summer' and 'winter'
## 
## Log intensity:  ~elev + I(elev^2) + marks
## 
## Fitted trend coefficients:
##   (Intercept)          elev     I(elev^2)   marksspring   markssummer 
## -1.850394e+01 -3.937116e-03  9.510112e-07 -3.957121e-02  7.891398e-01 
##   markswinter 
## -2.379296e+00 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -1.850394e+01 4.652449e-02 -1.859513e+01 -1.841276e+01   ***
## elev        -3.937116e-03 1.034097e-04 -4.139795e-03 -3.734437e-03   ***
## I(elev^2)    9.510112e-07 6.143974e-08  8.305915e-07  1.071431e-06   ***
## marksspring -3.957121e-02 5.053363e-02 -1.386153e-01  5.947288e-02      
## markssummer  7.891398e-01 4.266227e-02  7.055233e-01  8.727563e-01   ***
## markswinter -2.379296e+00 1.215116e-01 -2.617454e+00 -2.141137e+00   ***
##                     Zval
## (Intercept) -397.7247657
## elev         -38.0729907
## I(elev^2)     15.4787625
## marksspring   -0.7830669
## markssummer   18.4973701
## markswinter  -19.5808065

# Compares models using AIC
cat("AIC for poly model: ", AIC(model_poly_elevation), "\n")

## AIC for poly model:  142038.1

cat("AIC for gam model: ", AIC(gam_elevation), "\n")

## AIC for gam model:  4186.001

Model shows a negative relationship between bear locations and elevation, suggesting bears are less likely to be found at higher elevations. Quadratic term suggests a slight non-linear effect, but primary relationship appears negative. Model also indicates that summer sees the highest concentration of bear points, while winter has an expected decrease in bear points. Spring also shows a negative effect, although it is smaller than winter’s effect.

Bears tend to prefer lower elevations (based on the negative relationship with elevation). The summer season has the highest bear presence, while winter and spring have much lower bear occurrences, as to be expected.

Creates a point pattern object for each season

# Creates seasonal ppp objects
ppp_list <- list()
for (s in names(seasons)) {
  bears_loc_season <- bears_loc_filtered[bears_loc_filtered$month %in% seasons[[s]], ]
  bears_sf_season <- st_as_sf(bears_loc_season, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
  bears_sf_season_proj <- st_transform(bears_sf_season, crs = 3005)
  coords_season <- st_coordinates(bears_sf_season_proj)
  suppressWarnings(
    ppp_list[[s]] <- ppp(x = coords_season[,1], y = coords_season[,2], window = window)
  )
}

# Extracts covariates
elevation <- DATA$Elevation
cover <- DATA$Forest
hfi <- DATA$HFI
dist_water <- DATA$Dist_Water

Calculates and plots the Kernel Density Estimate for bear occurrences per each season

# KDE comparison
dens_list <- lapply(ppp_list, function(p) density(p, sigma = 10000))
for (s in names(dens_list)) {
  plot(dens_list[[s]], 
       main = paste("KDE Surface -", s))
}

Performs quadrat test for each season to analyze the spatial pattern of bear occurrences

# Quadrat test per season
for (s in names(ppp_list)) {
  qt <- quadrat.test(ppp_list[[s]], nx = 5, ny = 5)
  print(qt)
  plot(qt, 
       main = paste("Quadrat test -", s))
}

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 561.39, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 2202.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 3767, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp_list[[s]]
## X2 = 2720.5, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 21 tiles (irregular windows)

Dividing area into a grid of 5x5 and testing for deviations from a random distribution. Bear locations do not follow a random distribution and show spatial clustering across all seasons, meaning their presence is likely influenced by specific environmental or ecological factors that are season-dependent. The quadrant test results suggest that for all seasons (winter, spring, summer, and autumn), bear locations exhibit significant spatial clustering or non-random patterns rather than being randomly distributed (as indicated by the very small p-values).

Bear distribution is influenced by factors other than randomness, such as environmental features, resource availability, or behavior patterns that vary across seasons.

Winter: X2 = 561.39, p-value < 2.2e-16. P-value indicates a strong deviation from CSR for winter, suggesting that bear locations are not randomly distributed in the winter season.

Spring: X2 = 2202.7, p-value < 2.2e-16. Test statistic is very large, and the p-value is highly significant, indicating that the bear locations in spring also show non-random clustering or patterns.

Summer: X2 = 3767, p-value < 2.2e-16. Summer further supports a non-random pattern, with an even higher test statistic and significant p-value.

Autumn: X2 = 2720.5, p-value < 2.2e-16. Autumn shows a significant departure from CSR, confirming that bear locations are non-randomly distributed.

Fits an inhomogeneous Poisson ppm for each season using elevation. Model predicts intensity and plots for each season.

# Inhomogeneous Poisson model per season
for (s in names(ppp_list)) {
  model <- ppm(ppp_list[[s]] ~ elev, covariates = list(elev = elevation))
  plot(predict(model), 
       main = paste("Intensity ~ Elevation -", s))
}

## Warning: glm.fit: algorithm did not converge

## Warning: Values of the covariate 'elev' were NA or undefined at 0.02% (1 out of
## 5456) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev, data = NULL, interaction = NULL,

Fits a multivariate Poisson Point Process model to bear occurrences for each season, considering all covariates: elevation, forest cover, hfi, and distance to water. Predicts intensity and plots for each season.

# Multivariate ppm
for (s in names(ppp_list)) {
  model <- ppm(ppp_list[[s]] ~ elev + cover + hfi + dist_water,
               covariates = list(elev = elevation, cover = cover, hfi = hfi, dist_water = dist_water))
  plot(predict(model), 
       main = paste("Predicted Intensity -", s))
}

## Warning: Values of the covariate 'hfi' were NA or undefined at 0.68% (4 out of
## 587) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,

## Warning: glm.fit: algorithm did not converge

## Warning: Values of the covariate 'hfi' were NA or undefined at 0.45% (12 out of
## 2667) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,

## Warning: Values of the covariates 'elev', 'hfi' were NA or undefined at 0.4%
## (22 out of 5456) of the quadrature points. Occurred while executing: ppm.ppp(Q
## = ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,

## Warning: Values of the covariate 'hfi' were NA or undefined at 0.67% (18 out of
## 2698) of the quadrature points. Occurred while executing: ppm.ppp(Q =
## ppp_list[[s]], trend = ~elev + cover + hfi + dist_water,

For each season - assessing spatial distribution of bear occurrences (I.e., clustered, dispersed, or randomly distributed)

Ripley’s K with isotropic correction for statistical efficiency

# Ripley's K with Envelope for each season
for (s in names(ppp_list)) {
  K <- Kest(ppp_list[[s]])
  E <- envelope(ppp_list[[s]], 
                Kest, 
                correction = "border", 
                rank = 1, 
                nsim = 19, 
                fix.n = TRUE)
  
  plot(K, 
       main = paste("Ripley's K -", s))
  plot(E, 
       add = TRUE,
       lty = 2)  
}

## Generating 19 simulations of CSR with fixed number of points  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

## Generating 19 simulations of CSR with fixed number of points  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

## Generating 19 simulations of CSR with fixed number of points  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

## Generating 19 simulations of CSR with fixed number of points  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

Across all four seasons, the Ripley’s K function plots show a pattern of clustering. This can be concluded as the estimated K functions (iso, trans, bord) are generally significantly above the expected K function under complete spatial randomness and they often lie outside the simulation envelope.

Winter: Clustering may reflect limited movement due to snow, or congregation around scarce food sources during milder periods. It could also be indirectly related to dens.

Spring: Clustering could indicate post-hibernation dispersal patterns, concentration around newly available food sources, habitat preferences, or potential attraction to human-related food.

Summer: Clustering may reflect the distribution of abundant food resources, proximity to water sources, preferred habitats, and ongoing attraction to human-related food during peak bear activity.

Autumn: Clustering could reflect bears concentrating on late-season food sources for hibernation, movements related to finding den sites, and continued attraction to human-related food as natural food availability changes.